home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / rat3c.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  18.7 KB  |  587 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module rat3c)
  13.  
  14. ;;    THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 3.
  15. ;;    IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS
  16.  
  17. (DECLARE-TOP(GENPREFIX A_3))
  18.  
  19. (LOAD-MACSYMA-MACROS RATMAC)
  20.  
  21. (declare-top (special $float $keepfloat $algebraic $ratfac *alpha)
  22.      (special genvar))
  23.  
  24.  
  25. ;; List of GCD algorithms.  Default one is first.
  26. (DEFMVAR *GCDL* '($SUBRES $EZ $RED $SPMOD $MOD $SPHEN $EEZ $ALGEBRAIC))
  27.  
  28. (DEFMVAR $GCD (CAR *GCDL*))             ;Sparse Modular
  29.  
  30. (DEFUN CGCD (A B) (COND (MODULUS 1)
  31.             ((AND $KEEPFLOAT (OR (FLOATP A) (FLOATP B))) 1)
  32.             (T (GCD A B)))) 
  33.  
  34. (DEFMFUN PQUOTIENTCHK (A B) (IF (EQN B 1) A (PQUOTIENT A B)))
  35.  
  36. (DEFUN PTIMESCHK (A B) (COND ((EQN A 1) B) ((EQN B 1) A) (T (PTIMES A B)))) 
  37.  
  38. (DEFUN PFLOATP (X)
  39.        (CATCH 'FLOAT (COND ((PCOEFP X) (FLOATP X)) (T (PFLOATP1 X)))))
  40.  
  41. (DEFUN PFLOATP1 (X)
  42.   (MAPC #'(LAMBDA (Q) (COND ((PCOEFP Q)
  43.                  (COND ((FLOATP Q) (THROW 'FLOAT T))))
  44.                 ((PFLOATP1 Q))))
  45.     (CDR X))
  46.   NIL)
  47.  
  48. (DEFMFUN PGCD (X Y)
  49.        (SETQ X (CAR (PGCDA X Y NIL)))
  50.        (COND ((PMINUSP X) (PMINUS X))
  51.          (MODULUS (MONIZE X))
  52.          (T X)))
  53.  
  54. (DEFMFUN PLCM (X Y)
  55.        (SETQ X (PGCDCOFACTS X Y))
  56.        (PTIMES (CAR X) (PTIMES (CADR X) (CADDR X))))
  57.  
  58. (DEFUN PLCMCOFACTS (X Y)
  59.        (SETQ X (PGCDCOFACTS X Y))
  60.        (LIST (PTIMES (CAR X) (PTIMES (CADR X) (CADDR X)))
  61.          (CADDR X) (CADR X)))
  62.  
  63. (DEFUN PGCDCOFACTS (X Y) 
  64.        (LET ((A (PGCDA X Y T)))
  65.         (COND ((CDR A) A)
  66.           ((EQUAL (SETQ A (CAR A)) 1) (LIST 1 X Y))
  67.           ((AND $ALGEBRAIC (NOT (PCOEFP A)))
  68.            (CONS A (PROG2 (SETQ X (RQUOTIENT X A)
  69.                     Y (RQUOTIENT Y A)
  70.                     A (PGCDCOFACTS (CDR X) (CDR Y)))
  71.                   (LIST (PTIMES (CAR X) (CADDR A))
  72.                     (PTIMES (CAR Y) (CADR A))
  73.                     (PTIMES (CADR A) (CDR Y))))))
  74.           ((EQ A X) (LIST X 1 (PQUOTIENT Y X)))
  75.           ((EQ A Y) (LIST A (PQUOTIENT X Y) 1))
  76.           (T (LIST A (PQUOTIENT X A) (PQUOTIENT Y A))))))
  77.  
  78. (DEFUN PGCDA (X Y COFAC? &AUX A C) 
  79.   (COND ((NOT $GCD) (LIST 1 X Y))
  80.     ((AND $KEEPFLOAT (OR (PFLOATP X) (PFLOATP Y)))
  81.      (COND ((OR (PCOEFP X) (PCOEFP Y)
  82.             (PCOEFP (SETQ A (CAR (PTERMCONT X))))
  83.             (PCOEFP (SETQ A (PGCD A (CAR (PTERMCONT Y))))))
  84.         (LIST 1 X Y))
  85.            (T (LIST A))))
  86.     ((PCOEFP X)
  87.      (COND ((PCOEFP Y)
  88.         (CONS (SETQ A (CGCD X Y))
  89.               (AND COFAC?
  90.                (LIST (CQUOTIENT X A) ;(CQUOTIENT 0 0) = 0
  91.                  (CQUOTIENT Y A)))))
  92.            ((ZEROP X) (LIST Y X 1))
  93.            (T (LIST (PCONTENT1 (CDR Y) X)))))
  94.     ((PCOEFP Y) (COND ((ZEROP Y) (LIST X 1 Y))
  95.               (T (LIST (PCONTENT1 (CDR X) Y)))))
  96.     ((EQUAL X Y) (LIST X 1 1))
  97.     ($RATFAC (FPGCDCO X Y))
  98.     ((NOT (EQ (P-VAR X) (P-VAR Y)))
  99.      (LIST (IF (POINTERGP (P-VAR X) (P-VAR Y))
  100.            (OLDCONTENT1 (P-TERMS X) Y)
  101.            (OLDCONTENT1 (P-TERMS Y) X))))
  102.     ((PROGN (DESETQ (A X) (PTERMCONT X))
  103.         (DESETQ (C Y) (PTERMCONT Y))
  104.         (NOT (AND (EQUAL A 1) (EQUAL C 1))))
  105.      (MAPCAR #'PTIMES (MONOMGCDCO A C COFAC?) (PGCDA X Y COFAC?)))
  106.     ((AND (NOT $ALGEBRAIC) (NOT MODULUS)
  107.           (DESETQ (A . C) (LIN-VAR-FIND (NREVERSE (PDEGREEVECTOR X))
  108.                         (NREVERSE (PDEGREEVECTOR Y))
  109.                         (REVERSE GENVAR))))
  110.      (COND ((f= A 1) (LINHACK X Y (CAR C) (CADR C) COFAC?))
  111.            (T (SETQ A (LINHACK Y X A (CADR C) COFAC?))
  112.           (IF (CDR A) (RPLACD A (NREVERSE (CDR A))))
  113.           A)))
  114.     ((EQ $GCD '$SPMOD) (LIST (ZGCD X Y)))
  115.     ((EQ $GCD '$SUBRES) (LIST (OLDGCD X Y)))
  116.     ((EQ $GCD '$ALGEBRAIC)
  117.      (IF (OR (PALGP X) (PALGP Y))
  118.          (LET (($GCD '$SUBRES)) (LIST (OLDGCD X Y)))
  119.          (LET (($GCD '$SPMOD)) (LIST (ZGCD X Y)))))
  120.     ((EQ $GCD '$EZ) (EZGCD2 X Y))
  121.     ((EQ $GCD '$RED) (LIST (OLDGCD X Y)))
  122.     ((EQ $GCD '$MOD) (NEWGCD X Y MODULUS))
  123.     ((EQ $GCD '$SPHEN) (SPHGCD X Y))
  124.     ((EQ $GCD '$EEZ) (EEZGCD X Y))
  125.     ((NOT (MEMQ $GCD *GCDL*))
  126.      (MERROR "GCD set incorrectly:~%~M" $GCD))
  127.     (T (LIST 1 X Y))))
  128.  
  129. ;; (DEFUN PMINDEG (P) (IF (PCOEFP P) 0 (NXTTOLAST (CDR P))))
  130.  
  131. ;; (DEFUN PDEGRED (P N) (IF (ZEROP N) P (PQUOTIENT P (MAKE-POLY (P-VAR P) N 1))))
  132.  
  133. (DEFUN MONOMGCDCO (P Q COFAC?)
  134.   (LET ((GCD (MONOMGCD P Q)))
  135.     (CONS GCD (IF COFAC? (LIST (PQUOTIENT P GCD) (PQUOTIENT Q GCD)) ()))))
  136.   
  137. (DEFUN MONOMGCD (P Q)
  138.   (COND ((OR (PCOEFP P) (PCOEFP Q)) 1)
  139.     ((EQ (P-VAR P) (P-VAR Q))
  140.      (MAKE-POLY (P-VAR P) (MIN (P-LE P) (P-LE Q))
  141.             (MONOMGCD (P-LC P) (P-LC Q))))
  142.     ((POINTERGP (CAR P) (CAR Q)) (MONOMGCD (P-LC P) Q))
  143.     (T (MONOMGCD P (P-LC Q)))))
  144.     
  145.  
  146. (DEFUN LINHACK (POL1 POL2 NONLINDEG VAR COFAC?)
  147.   (PROG (COEFF11 COEFF12 GCDAB RPOL1 RPOL2 GCDCD GCDCOEF)
  148.     (DESETQ (COEFF11 . COEFF12) (BOTHPRODCOEF (MAKE-POLY VAR) POL1))
  149.     (SETQ GCDAB (IF (PZEROP COEFF12) COEFF11
  150.             (PGCD COEFF11 COEFF12)))
  151.     (COND ((EQUAL GCDAB 1)
  152.            (COND ((SETQ COEFF11 (TESTDIVIDE POL2 POL1))
  153.               (RETURN (LIST POL1 1 COEFF11)))
  154.              (T (RETURN (LIST 1 POL1 POL2))))))
  155.     (SETQ RPOL1 (PQUOTIENT POL1 GCDAB))
  156.     (DESETQ (GCDCD RPOL2) (LINHACKCONTENT VAR POL2 NONLINDEG))
  157.     (COND ((EQUAL GCDCD 1)
  158.            (COND ((SETQ COEFF12 (TESTDIVIDE RPOL2 RPOL1))
  159.               (RETURN (LIST RPOL1 GCDAB COEFF12)))
  160.              (T (RETURN (LIST 1 POL1 POL2))))))
  161.     (COND (COFAC? (DESETQ (GCDCOEF COEFF11 COEFF12)
  162.                   (PGCDCOFACTS GCDAB GCDCD))
  163.               (COND ((SETQ GCDCD (TESTDIVIDE RPOL2 RPOL1))
  164.                  (RETURN (LIST (PTIMES GCDCOEF RPOL1)
  165.                        COEFF11
  166.                        (PTIMES COEFF12 GCDCD))))
  167.                 (T (RETURN (LIST GCDCOEF 
  168.                          (PTIMES COEFF11 RPOL1)
  169.                          (PTIMES COEFF12 RPOL2))))))
  170.           (T (SETQ GCDCOEF (PGCD GCDCD GCDAB))
  171.          (COND ((TESTDIVIDE RPOL2 RPOL1)
  172.             (RETURN (LIST (PTIMES GCDCOEF RPOL1))))
  173.                (T (RETURN (LIST GCDCOEF))))))))
  174.  
  175. (DEFUN LIN-VAR-FIND (A B C)
  176.   (DO ((VARL C (CDR VARL))
  177.        (DEGL1 A (CDR DEGL1))
  178.        (DEGL2 B (CDR DEGL2)))
  179.       ((OR (NULL DEGL1) (NULL DEGL2)) NIL)
  180.     (IF (EQUAL (MIN (CAR DEGL1) (CAR DEGL2)) 1)
  181.     (RETURN (LIST (CAR DEGL1) (CAR DEGL2) (CAR VARL))))))
  182.  
  183. (DEFUN LINHACKCONTENT (VAR POL NONLINDEG &AUX (NPOL POL) COEF GCD)
  184.   (DO ((I NONLINDEG (f1- I)))
  185.       ((f= I 0) (LIST (SETQ GCD (PGCD GCD NPOL)) (PQUOTIENT POL GCD)))
  186.     (DESETQ (COEF . NPOL) (BOTHPRODCOEF (MAKE-POLY VAR I 1) NPOL)) 
  187.     (UNLESS (PZEROP COEF)
  188.         (SETQ GCD (IF (NULL GCD) COEF (PGCD COEF GCD)))
  189.         (IF (EQUAL GCD 1) (RETURN (LIST 1 POL))))))
  190.  
  191. ;;*** THIS IS THE REDUCED POLYNOMIAL REMAINDER SEQUENCE GCD (COLLINS')
  192.  
  193. (DEFUN OLDGCD (X Y &AUX U V S EGCD)            ;only called from pgcda
  194.   (DESETQ (X  U) (OLDCONTENT X))
  195.   (DESETQ (Y  V) (OLDCONTENT Y))
  196.   (SETQ EGCD (GCD (PGCDEXPON U) (PGCDEXPON V)))
  197.   (IF (> EGCD 1)
  198.       (SETQ U (PEXPON*// U EGCD NIL)
  199.         V (PEXPON*// V EGCD NIL)))
  200.   (IF (f> (P-LE V) (P-LE U)) (EXCH U V))
  201.   (SETQ S (CASE $GCD
  202.         ($RED (REDGCD U V))
  203.         ($SUBRES (SUBRESGCD U V))
  204.         (T (MERROR "Illegal GCD algorithm"))))
  205.   (UNLESS (EQUAL S 1)
  206.       (SETQ S (PEXPON*// (PRIMPART
  207.                   (IF $ALGEBRAIC S
  208.                   (PQUOTIENT S (PQUOTIENT (P-LC S)
  209.                               (PGCD (P-LC U)
  210.                                 (P-LC V))))))
  211.                  EGCD T)))
  212.   (SETQ S (PTIMESCHK S (PGCD X Y)))
  213.   (AND $ALGEBRAIC (NOT (PCOEFP (SETQ U (LEADALGCOEF S))))
  214.        (NOT (EQUAL U S)) (SETQ S (ALGNORMAL S)))
  215.   (COND (MODULUS (MONIZE S))
  216.     ((PMINUSP S) (PMINUS S)) 
  217.     (T S)))
  218.  
  219. (DEFUN PGCDEXPON (P)
  220.   (IF (PCOEFP P) 0
  221.       (DO ((D (CADR P) (GCD D (CAR L)))
  222.        (L (CDDDR P) (CDDR L)))
  223.       ((OR (NULL L) (f= D 1)) D))))
  224.  
  225. (DEFUN PEXPON*// (P N *?)
  226.   (IF (OR (PCOEFP P) (f= N 1)) P
  227.       (DO ((ANS (LIST (CAR P))
  228.         (CONS (CADR L)
  229.               (CONS (IF *? (f* (CAR L) N)
  230.                 (// (CAR L) N))
  231.                 ANS)))
  232.        (L (CDR P) (CDDR L)))
  233.       ((NULL L) (NREVERSE ANS)))))
  234.  
  235. ;;polynomial gcd using reduced prs
  236.  
  237. (defun redgcd (p q &aux (d 0))
  238.   (sloop until (zerop (pdegree q (p-var p)))
  239.     do (psetq p q
  240.           q (pquotientchk (prem p q) (pexpt (p-lc p) d))
  241.           d (f+ (p-le p) 1 (f- (p-le q))))
  242.     finally (return (if (pzerop q) p 1))))
  243.  
  244. ;;computes gcd's using subresultant prs TOMS Sept. 1978
  245.  
  246. (defun subresgcd (p q)                    
  247.   (sloop for g = 1 then (p-lc p)
  248.     for h = 1 then (pquotient (pexpt g d) h^1-d)
  249.     for d = (f- (p-le p) (p-le q))
  250.     for h^1-d = (if (equal h 1) 1 (pexpt h (f1- d)))
  251.         do (psetq p q
  252.           q (pquotientchk (prem p q) (ptimes g (ptimes h h^1-D))))
  253.     if (zerop (pdegree q (p-var p))) return (if (pzerop q) p 1)))
  254.  
  255. ;;*** THIS COMPUTES PSEUDO REMAINDERS
  256.  
  257. ;(DECLARE (SPECIAL K LCU LCV) (FIXNUM K M I))
  258.  
  259. (DEFUN PSQUOREM1 (U V QUOP)
  260.        (PROG (K (M 0) LCU LCV QUO LC)
  261.          (declare (special k lcu lcv) (fixnum #-cl k m))
  262.          (SETQ LCV (PT-LC V))
  263.          (SETQ K (f- (PT-LE U) (PT-LE V)))
  264.          (COND ((MINUSP K) (RETURN (LIST 1 '(0 0) U))))
  265.          (IF QUOP (SETQ LC (PEXPT (PT-LC V) (f1+ K))))
  266.        A     (SETQ LCU (PMINUS (PT-LC U)))
  267.              (IF QUOP (SETQ QUO (CONS (PTIMES (PT-LC U) (PEXPT (PT-LC V) K))
  268.                       (CONS K QUO))))
  269.          (COND ((NULL (SETQ U (PGCD2 (PT-RED U) (PT-RED V))))
  270.             (RETURN (LIST LC (NREVERSE QUO) '(0 0))))
  271.            ((MINUSP (SETQ M (f- (PT-LE U) (PT-LE V))))
  272.             (SETQ U (COND ((ZEROP K) U)
  273.                   (T (PCTIMES1 (PEXPT LCV K) U))))
  274.             (RETURN (LIST LC (NREVERSE QUO) U)))
  275.            ((f> (f1- K) M)
  276.             (SETQ U (PCTIMES1 (PEXPT LCV (f- (f1- K) M)) U))))
  277.          (SETQ K M)
  278.          (GO A)))
  279.  
  280. (DEFUN PREM (P Q)
  281.   (COND ((PCOEFP P) (IF (PCOEFP Q) (CREMAINDER P Q) P))
  282.     ((PCOEFP Q) (PZERO))
  283.     (T (PSIMP (P-VAR P) (PGCD1 (P-TERMS P) (P-TERMS Q))))))
  284.  
  285. (DEFMFUN PGCD1 (U V) (CADDR (PSQUOREM1 U V NIL)))
  286.  
  287. (DEFUN PGCD2 (U V &AUX (I 0))
  288.   (declare (special k lcu lcv) (fixnum #-cl k i))
  289.   (COND ((NULL U) (PCETIMES1 V K LCU))
  290.     ((NULL V) (PCTIMES1 LCV U))
  291.     ((ZEROP (SETQ I (f+ (PT-LE U) (f- K) (f- (CAR V)))))
  292.      (PCOEFADD (PT-LE U) (PPLUS (PTIMES LCV (PT-LC U))
  293.                     (PTIMES LCU (PT-LC V)))
  294.            (PGCD2 (PT-RED U) (PT-RED V))))
  295.     ((MINUSP I)
  296.      (LIST* (f+ (PT-LE V) K) (PTIMES LCU (PT-LC V)) (PGCD2 U (PT-RED V))))
  297.     (T (LIST* (PT-LE U) (PTIMES LCV (PT-LC U)) (PGCD2 (PT-RED U) V)))))
  298.  
  299. ;(DECLARE (UNSPECIAL K LCU LCV) (NOTYPE K M I))
  300.        
  301. ;;;*** OLDCONTENT REMOVES ALL BUT MAIN VARIABLE AND PUTS THAT IN CONTENT
  302. ;;;***  OLDCONTENT OF 3*A*X IS 3*A (WITH MAINVAR=X)
  303.  
  304. (DEFUN RCONTENT (P)    ;RETURNS RAT-FORMS
  305.        (LET ((Q (OLDCONTENTA P)))
  306.         (LIST (CONS Q 1) (COND ($ALGEBRAIC (RQUOTIENT P Q))
  307.                    (T (CONS (PQUOTIENT P Q) 1))))))
  308.  
  309. (DEFUN OLDCONTENTA (X)
  310.        (COND ((PCOEFP X) X)
  311.          (T (SETQ X (CONTSORT (CDR X)))
  312.         (OLDCONTENT2 (CDR X) (CAR X)))))
  313.  
  314. (DEFMFUN OLDCONTENT (X)
  315.        (COND ((PCOEFP X) (LIST X 1))
  316.          ((NULL (P-RED X))
  317.           (LIST (P-LC X) (MAKE-POLY (P-VAR X) (P-LE X) 1)))
  318.          (T (LET ((U (CONTSORT (CDR X))) V)
  319.                 (SETQ U (OLDCONTENT2 (CDR U) (CAR U))
  320.                       V (COND ($ALGEBRAIC (CAR (RQUOTIENT X U)))
  321.                   (T (PCQUOTIENT X U))))
  322.                 (COND ((PMINUSP V) (LIST (PMINUS U) (PMINUS V)))
  323.                       (T (LIST U V)))))))
  324.  
  325. (DEFUN OLDCONTENT1 (X GCD) (COND ((EQUAL GCD 1) 1)
  326.                  ((NULL X) GCD)
  327.                  (T (OLDCONTENT2 (CONTSORT X) GCD))))
  328.  
  329. (DEFUN OLDCONTENT2 (X GCD)
  330.        (do ((x x (cdr x))
  331.         (gcd gcd (pgcd (car x) gcd)))
  332.        ((or (null x) (equal gcd 1)) gcd)))
  333.  
  334. (DEFUN CONTSORT (X)
  335.        (SETQ X (COEFL X))
  336.        (COND ((zl-MEMBER 1 X)'(1))
  337.          ((NULL (CDR X))X)
  338.          (T (SORT X (FUNCTION CONTODR)))))
  339.  
  340. (DEFUN COEFL (X)
  341.        (do ((x x (cddr x))
  342.         (ans nil (cons (cadr x) ans)))
  343.        ((null x) ans)))
  344.  
  345. (DEFUN CONTODR (A B)
  346.        (COND ((PCOEFP A) T)
  347.          ((PCOEFP B) NIL)
  348.          ((EQ (CAR A) (CAR B)) (NOT (f> (CADR A) (CADR B))))
  349.          (T (POINTERGP (CAR B)(CAR A)))))
  350.  
  351. ;;;*** PCONTENT COMPUTES INTEGER CONTENT
  352. ;;;*** PCONTENT OF 3*A*X IS 3 IF MODULUS = NIL  1 OTHERWISE
  353.  
  354. (DEFUN PCONTENT (X)
  355.        (COND ((PCOEFP X) (LIST X 1))
  356.          (T (LET ((U (PCONTENTZ X)))
  357.           (IF (EQN U 1) (LIST 1 X)
  358.               (LIST U (PCQUOTIENT X U)))))))
  359.  
  360. (DEFUN PCONTENT1 (X GCD)
  361.        (DO ((X X (CDDR X))
  362.         (GCD GCD (CGCD GCD (PCONTENTZ (CADR X)))))
  363.        ((OR (NULL X) (EQUAL GCD 1)) GCD)))
  364.  
  365. (DEFUN PCONTENTZ (P)
  366.        (COND ((PCOEFP P) P)
  367.          (T (PCONTENT1 (P-RED P) (PCONTENTZ (P-LC P))))))
  368.  
  369. (DEFUN UCONTENT (P)                    ;CONTENT OF UNIV. POLY
  370.   (COND ((PCOEFP P) (ABS P))
  371.     (T (SETQ P (MAPCAR #'ABS (COEFL (CDR P))))
  372.        (LET ((M (APPLY #'MIN P)))
  373.         (OLDCONTENT2 (zl-DELETE M P) M)))))
  374.  
  375. ;;***    PGCDU CORRESPONDS TO BROWN'S ALGORITHM U
  376.  
  377. ;;;PGCDU IS NOT NOW IN RAT;UFACT >
  378.  
  379. (DEFMFUN PGCDU (P Q)
  380.   (DO () ((PZEROP Q) (MONIZE P))
  381.     (PSETQ P Q Q (PMODREM P Q))))
  382.  
  383. ;(DECLARE (SPECIAL K Q* QUO) (FIXNUM K))
  384.  
  385. (DEFUN PMODREM (X Y)
  386.   (COND ((NULL MODULUS)
  387.      (MERROR "Illegal use of PMODREM"))
  388.     ((PACOEFP Y) (IF (PZEROP Y) X 0))
  389.     ((PACOEFP X) X)
  390.     ((EQ (P-VAR X) (P-VAR Y))
  391.      (PSIMP (CAR X) (PGCDU1 (P-TERMS X) (P-TERMS Y) NIL)))
  392.     (T (MERROR "Illegal use of PMODREM"))))
  393.  
  394. (DEFUN PMODQUO (U V &AUX QUO)
  395.   (declare (special quo))
  396.   (COND ((NULL MODULUS)
  397.      (MERROR "Illegal use of PMODQUO"))
  398.     ((PCOEFP V) (CONS (PTIMES (CRECIP V) U) 0))
  399.     ((ALG V) (CONS (PTIMES (PAINVMOD V) U) 0))
  400.     ((PACOEFP U) (CONS 0 U))
  401.     ((NOT (EQ (P-VAR U) (P-VAR V)))
  402.      (MERROR "Illegal use of PMODQUO"))
  403.     (T (XCONS (PSIMP (CAR U) (PGCDU1 (CDR U) (CDR V) T))
  404.           (PSIMP (CAR U) QUO)))))
  405.  
  406. (COMMENT (DEFUN PMODREM (X Y) (COND ((NULL MODULUS)
  407.                   (MERROR "Illegal use of PMODREM"))
  408.                  ((PACOEFP Y) 0)
  409.                  ((PACOEFP X) X)
  410.                  ((EQ (CAR X) (CAR Y))
  411.                   (DPDISREP
  412.                    (DPMODREM (DPREP X) (DPREP Y))))
  413.                  (T (MERROR "Illegal use of PMODREM")))))
  414.  
  415. (DEFUN PGCDU1 (U V PQUO*)
  416.   (let ((invv (PAINVMOD (PT-LC V))) (k 0) q*)
  417.     (declare (special k quo q*) (fixnum k))
  418.     (SLOOP UNTIL (MINUSP (setq K (f- (PT-LE U) (PT-LE V))))
  419.       DO (SETQ Q* (PTIMES INVV (PT-LC U)))
  420.       IF PQUO* DO (SETQ QUO (NCONC QUO (LIST K Q*)))
  421.       WHEN (PTZEROP (SETQ U (PQUOTIENT2 (PT-RED U) (PT-RED V))))
  422.        RETURN (PTZERO)
  423.       FINALLY (RETURN U))))
  424.  
  425. ;(DECLARE (UNSPECIAL K Q* QUO) (NOTYPE K))
  426.  
  427.  
  428. (DEFUN NEWPRIME (P)
  429.   (declare (special bigprimes))        ;defined later on
  430.   (COND ((NULL P) (CAR BIGPRIMES))
  431.     (T (DO ((PL BIGPRIMES (CDR PL)))
  432.            ((NULL PL) (SETQ P (FNEWPRIME P))
  433.               (SETQ BIGPRIMES (NCONC BIGPRIMES (LIST P)))
  434.               P)
  435.          (IF (f< (CAR PL) P) (RETURN (CAR PL)))))))
  436.  
  437. (DEFUN FNEWPRIME (P)    ; Finds biggest prime less than fixnum P.
  438.   (DO ((PP (IF (ODDP P) (f- P 2) (f- P 1)) (f- PP 2))) ((f< PP 0))
  439.       (IF (PRIMEP PP) (RETURN PP))))
  440.  
  441. (DEFUN PRIMEP (P)
  442.   (AND (OR (LESSP P 14.)
  443.        (LET ((MODULUS P))
  444.         (AND (EQUAL 1 (CEXPT 13. (SUB1 P))) (EQUAL 1 (CEXPT 3 (SUB1 P))))))
  445.        (NULL (CDDR (SETQ P (CFACTORW P))))
  446.        (= 1 (CADR P))))
  447.  
  448. ;; #O <form> reads <form> in octal (base 8)
  449.  
  450.  
  451. (DEFVAR BIGPRIMES NIL)
  452. #+CL
  453. (eval-when (load )
  454.  
  455. ;; it is convenient to have the bigprimes be actually less than
  456. ;; half the size of the most positive fixnum, so that arithmetic is
  457. ;; easier
  458. #.
  459. (case most-positive-fixnum
  460.   (2147483647
  461.     '(setq bigprimes
  462.        '(1073741789 1073741783 1073741741 1073741723 1073741719 1073741717
  463.  1073741689 1073741671 1073741663 1073741651 1073741621 1073741567
  464.  1073741561 1073741527 1073741503 1073741477 1073741467 1073741441
  465.  1073741419 1073741399)
  466. ))
  467.   ;; Could always use the following, but it takes several seconds to compute
  468.   ;; so if we want to autoload this file, it is tiresome.
  469.   (t '(DO ((I 0 (f1+ I))                ;GENERATES 20 LARGEST
  470.      (P (quotient most-positive-fixnum 2) (NEWPRIME P)))        ;PRIMES < WORD
  471.     ((= I 20.)))))
  472.  
  473. (setq *ALPHA (CAR BIGPRIMES))
  474. )
  475.  
  476. (DEFMVAR *ALPHA (CAR BIGPRIMES))
  477.  
  478. #+MacLisp 
  479. (DEFMVAR BIGPRIMES
  480.       '(#O 377777777741 #O 377777777717 #O 377777777703 #O 377777777673
  481.     #O 377777777661 #O 377777777607 #O 377777777563 #O 377777777411
  482.     #O 377777777313 #O 377777777273 #O 377777777233 #O 377777777075
  483.     #O 377777777015 #O 377777776771 #O 377777776755 #O 377777776735
  484.     #O 377777776725 #O 377777776677 #O 377777776661 #O 377777776653))
  485.  
  486.  
  487. ;;; list of primes less than 2^30 -1  (limit of smallnums in Franz/vax)
  488. #+Franz 
  489. (defmvar bigprimes '(1073741789. 1073741783. 1073741741. 1073741723. 
  490.           1073741719. 1073741717. 1073741689. 1073741671.
  491.           1073741663. 1073741651. 1073741621. 1073741567. 
  492.           1073741561. 1073741527. 1073741503. 1073741477.
  493.           1073741467. 1073741441. 1073741419. 1073741399.))
  494.  
  495. #+NIL
  496. (PROGN 'COMPILE
  497. ;;; It takes a lot longer to compute the 35 bit primes for the PDP-10,
  498. ;;; thats why they are wired-in there.
  499. (DEFMVAR BIGPRIMES
  500.   '(#o3777777775 #o3777777737 #o3777777725 #o3777777701 #o3777777667
  501.     #o3777777665 #o3777777643 #o3777777635 #o3777777607 #o3777777573
  502.     #o3777777557 #o3777777527 #o3777777511 #o3777777503 #o3777777475
  503.     #o3777777455 #o3777777433 #o3777777401 #o3777777361 #o3777777343))
  504.  
  505. (DEFUN BIGPRIMES ()
  506.   (SETQ BIGPRIMES ())
  507.   (DO ((I 0 (f1+ I)) (P MOST-POSITIVE-FIXNUM (NEWPRIME P)))
  508.       ((= I 20.))))
  509. ;; wire it in for now. PRIMEP is losing. (BIGPRIMES)
  510. )
  511.  
  512. ;#+CL
  513. ;(DO ((I 0 (f1+ I))                ;GENERATES 20 LARGEST
  514. ;         (P (LSH -1 -1) (NEWPRIME P)))        ;PRIMES < WORD
  515. ;        ((= I 20.)))
  516.  
  517.  
  518.  
  519. (DEFMFUN $PRIMEP (P)
  520.  (IF (NOT (INTEGERP P)) (MERROR "Argument to PRIMEP must be an integer:~%~M" P))
  521.  (LET ($INTFACLIM) (PRIMEP (ABS P))))
  522.  
  523.  
  524.  
  525. (DEFUN LEADCOEFFICIENT (P) (IF (PCOEFP P) P (LEADCOEFFICIENT (CADDR P))))
  526.  
  527. (DEFUN MAXCOEFFICIENT (P) (IF (PCOEFP P) (ABS P) (MAXCOEF1 (CDR P))))
  528.  
  529. (DEFUN MAXCOEF1 (P)
  530.   (IF (NULL P) 0 (MAX (MAXCOEFFICIENT (CADR P)) (MAXCOEF1 (CDDR P)))))
  531.  
  532. (DEFUN MAXNORM (POLY)
  533.   (IF (NULL POLY) 0 (MAX (NORM (CADR POLY)) (MAXNORM (CDDR POLY)))))
  534.  
  535. (DEFUN NORM (POLY)
  536.        (COND ((NULL POLY) 0)
  537.          ((PCOEFP POLY) (ABS POLY))
  538.          (T (PLUS (NORM (CADDR POLY)) (NORM1 (CDDDR POLY)) )) ))
  539.  
  540. (DEFUN NORM1 (POLY)
  541.   (IF (NULL POLY) 0 (PLUS (NORM (CADR POLY)) (NORM1 (CDDR POLY)) )) )
  542.  
  543. (DEFMFUN PDEGREE (P VAR)
  544.   (COND ((PCOEFP P) 0)
  545.     ((EQ VAR (P-VAR P)) (P-LE P))
  546.     ((POINTERGP VAR (P-VAR P)) 0)
  547.     (T (DO ((L (P-RED P) (PT-RED L))
  548.         (E (PDEGREE (P-LC P) VAR) (MAX E (PDEGREE (PT-LC L) VAR))))
  549.            ((NULL L) E)))))
  550.  
  551. (DEFUN POLY-IN-VAR (P V)
  552.   (COND ((OR (PCOEFP P) (POINTERGP V (P-VAR P))) (LIST 0 P))
  553.     ((EQ (P-VAR P) V) (P-TERMS P))
  554.     ((SLOOP WITH ANS
  555.            FOR (EXP COEF) ON (P-TERMS P) BY 'PT-RED
  556.            DO (SETQ ANS (PPLUS1 ANS
  557.                     (EVERYSUBST2 (POLY-IN-VAR COEF V)
  558.                          (LIST (P-VAR P) EXP 1))))
  559.            FINALLY (RETURN ANS)))))
  560.  
  561. (DEFUN UNIVAR (X) (OR (NULL X) (AND (PCOEFP (PT-LC X)) (UNIVAR (PT-RED X)))))
  562.  
  563. ;;**THE CHINESE REMAINDER ALGORITHM IS A SPECIAL CASE OF LAGRANGE INTERPOLATION
  564.  
  565. (DEFUN LAGRANGE3 (U UK P QK)
  566.        (SETQMODULUS P)
  567.        (SETQ UK (PDIFFERENCE UK (PMOD U)))
  568.        (COND ((PZEROP UK) (SETQ MODULUS NIL) U)
  569.          (T (SETQ UK (PCTIMES (CRECIP (CMOD QK)) UK))
  570.         (SETQ MODULUS NIL)
  571.         (PPLUS U (PCTIMES QK UK)))))
  572.  
  573.  
  574. (DEFUN LAGRANGE33 (U UK QK XK)
  575.   (DECLARE (SPECIAL XV))
  576.   (SETQ UK (PDIFFERENCE UK (PCSUBST U XK XV)))
  577.   (COND ((PZEROP UK) U)
  578.     (T (PPLUS U (PTIMES
  579.              (PCTIMES (CRECIP (PCSUBST QK XK XV)) UK)
  580.              QK)))))
  581.  
  582.  
  583. ;;;*************************************************************
  584.  
  585. ;;    THIS IS THE END THE NEW RATIONAL FUNCTION PACKAGE PART 3.
  586. ;;    IT INCLUDES THE GCD ROUTINES, THEIR SUPPORTING FUNCTIONS
  587.